# clusteringFunUtils.R
# 
# Author: yaping
# 2013
# Aug 6, 2013
# 10:33:33 AM
###############################################################################



##default function to use:
##in nome-seq mode and two-step clustering mode, this input files' first 4 should be HCT116 methy, DKO1 methy, HCT116 acc, DKO1 acc 
##reorder the matrix by both of samples' DNA methylation, then by DKO's accessibility level
clusteringNOMeSeqMatrix1step<-function(methy_HCT116,methy_DKO1, accSubClusterNum=4){
	content_hcg_HCT116<-read.table(methy_HCT116,sep="\t",header=F)	
	content_hcg_HCT116<-content_hcg_HCT116[!duplicated(content_hcg_HCT116),]
	rownames(content_hcg_HCT116)<-paste(content_hcg_HCT116[,1],content_hcg_HCT116[,2],content_hcg_HCT116[,3],content_hcg_HCT116[,4],sep=":")
	
	content_hcg_DKO1<-read.table(methy_DKO1,sep="\t",header=F)	
	content_hcg_DKO1<-content_hcg_DKO1[!duplicated(content_hcg_DKO1),]
	rownames(content_hcg_DKO1)<-paste(content_hcg_DKO1[,1],content_hcg_DKO1[,2],content_hcg_DKO1[,3],content_hcg_DKO1[,4],sep=":")
	
	
	dataSeq<-seq(5+((length(content_hcg_HCT116[1,])-5)/2)-as.integer(scale/bin_size_align), 5+((length(content_hcg_HCT116[1,])-5)/2)+as.integer(scale/bin_size_align), by=as.integer(move_step/bin_size_align))
	value_hcg_HCT116<-NULL
	value_hcg_DKO1<-NULL
	
	for(i in dataSeq){
		
		if(floor(i+move_step/(2*bin_size_align)-1)-ceiling(i-move_step/(2*bin_size_align)) <= 0){
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align))]))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
		}else{
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(rowMeans(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(rowMeans(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))
		}
		
	}
	
	rownames(value_hcg_HCT116)<-rownames(content_hcg_HCT116)
	rownames(value_hcg_DKO1)<-rownames(content_hcg_DKO1)
	
	
	##clustering:
	valueTmpMethyNoCluster.1<-value_hcg_DKO1[,regionToCluster]
	valueTmpMethyNoCluster.1<-interpNA(valueTmpMethyNoCluster.1[rowSums(!is.na(valueTmpMethyNoCluster.1))>=3,], method = "linear")
	valueTmpMethyNoCluster.1<-removeNA(valueTmpMethyNoCluster.1)
	
	valueTmpMethyNoCluster.2<-value_hcg_HCT116[,regionToCluster]
	valueTmpMethyNoCluster.2<-interpNA(valueTmpMethyNoCluster.2[rowSums(!is.na(valueTmpMethyNoCluster.2))>=3,], method = "linear")
	valueTmpMethyNoCluster.2<-removeNA(valueTmpMethyNoCluster.2)
	common.methy<-intersect(rownames(valueTmpMethyNoCluster.1),rownames(valueTmpMethyNoCluster.2))
	valueTmpMethyNoCluster<-cbind(valueTmpMethyNoCluster.1[common.methy,],valueTmpMethyNoCluster.2[common.methy,])
	
	#valueTmpAccNoCluster<-cbind(value_hcg_DKO1[,regionToCluster],value_hcg_HCT116[,regionToCluster])
	#valueTmpAccNoCluster<-valueTmpAccNoCluster[rowSums(is.na(valueTmpAccNoCluster))<=quantile(rowSums(is.na(valueTmpAccNoCluster)),probs=seq(0,1,by=0.05))[19],]
	
	cluster.acc<-hclust(dist(valueTmpMethyNoCluster),method="ward")
	subClusterOrder<-cutree(cluster.acc, k = accSubClusterNum)
	names(subClusterOrder)<-rownames(valueTmpMethyNoCluster)
	combinedRowOrder<-rownames(valueTmpMethyNoCluster[cluster.acc$order,])
	
	
	returnResult<-list(rowNames=combinedRowOrder, subClusterOrder=subClusterOrder, capUplimit=NULL, capDownLimit=NULL, dendgram=as.dendrogram(cluster.acc))
	return(returnResult)
	
}

clusteringNOMeSeqMatrix1step<-function(methy_DKO1, accSubClusterNum=4){
	
	content_hcg_DKO1<-read.table(methy_DKO1,sep="\t",header=F)	
	content_hcg_DKO1<-content_hcg_DKO1[!duplicated(content_hcg_DKO1),]
	rownames(content_hcg_DKO1)<-paste(content_hcg_DKO1[,1],content_hcg_DKO1[,2],content_hcg_DKO1[,3],content_hcg_DKO1[,4],sep=":")
	
	
	dataSeq<-seq(5+((length(content_hcg_DKO1[1,])-5)/2)-as.integer(scale/bin_size_align), 5+((length(content_hcg_DKO1[1,])-5)/2)+as.integer(scale/bin_size_align), by=as.integer(move_step/bin_size_align))
	value_hcg_DKO1<-NULL
	
	for(i in dataSeq){
		
		if(floor(i+move_step/(2*bin_size_align)-1)-ceiling(i-move_step/(2*bin_size_align)) <= 0){
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
		}else{
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(rowMeans(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))
		}
		
	}
	
	rownames(value_hcg_DKO1)<-rownames(content_hcg_DKO1)
	
	
	##clustering:
	#valueTmpMethyNoCluster.1<-value_hcg_DKO1[,regionToCluster]
	#valueTmpMethyNoCluster.1<-interpNA(valueTmpMethyNoCluster.1[rowSums(!is.na(valueTmpMethyNoCluster.1))>=3,], method = "linear")
	#valueTmpMethyNoCluster.1<-removeNA(valueTmpMethyNoCluster.1)
	#valueTmpMethyNoCluster<-valueTmpMethyNoCluster.1
	valueTmpMethyNoCluster<-value_hcg_DKO1[,regionToCluster]	
	#maxNAallow=quantile(rowSums(is.na(valueTmpMethyNoCluster)),probs=seq(0,1,0.05),na.rm =T)[19]
	maxNAallow=as.integer(length(regionToCluster)/2)
	#maxNAallow=5
	valueTmpMethyNoCluster<-valueTmpMethyNoCluster[rowSums(is.na(valueTmpMethyNoCluster))<=maxNAallow,]
	
	
	
	
	#valueTmpAccNoCluster<-cbind(value_hcg_DKO1[,regionToCluster],value_hcg_HCT116[,regionToCluster])
	#valueTmpAccNoCluster<-valueTmpAccNoCluster[rowSums(is.na(valueTmpAccNoCluster))<=quantile(rowSums(is.na(valueTmpAccNoCluster)),probs=seq(0,1,by=0.05))[19],]
	
	cluster.acc<-hclust(dist(valueTmpMethyNoCluster),method="ward")
	subClusterOrder<-cutree(cluster.acc, k = accSubClusterNum)
	names(subClusterOrder)<-rownames(valueTmpMethyNoCluster)
	combinedRowOrder<-rownames(valueTmpMethyNoCluster[cluster.acc$order,])
	
	
	returnResult<-list(rowNames=combinedRowOrder, subClusterOrder=subClusterOrder, capUplimit=NULL, capDownLimit=NULL, dendgram=as.dendrogram(cluster.acc))
	return(returnResult)
	
}

clusteringNOMeSeqMatrix2step<-function(methy_HCT116,methy_DKO1,acc_HCT116, acc_DKO1, accSubClusterNum=4, methySubClusterNum=2){
	content_hcg_HCT116<-read.table(methy_HCT116,sep="\t",header=F)	
	content_hcg_HCT116<-content_hcg_HCT116[!duplicated(content_hcg_HCT116),]
	rownames(content_hcg_HCT116)<-paste(content_hcg_HCT116[,1],content_hcg_HCT116[,2],content_hcg_HCT116[,3],content_hcg_HCT116[,4],sep=":")
	
	content_hcg_DKO1<-read.table(methy_DKO1,sep="\t",header=F)	
	content_hcg_DKO1<-content_hcg_DKO1[!duplicated(content_hcg_DKO1),]
	rownames(content_hcg_DKO1)<-paste(content_hcg_DKO1[,1],content_hcg_DKO1[,2],content_hcg_DKO1[,3],content_hcg_DKO1[,4],sep=":")
	
	content_gch_HCT116<-read.table(acc_HCT116,sep="\t",header=F)	
	content_gch_HCT116<-content_gch_HCT116[!duplicated(content_gch_HCT116),]
	rownames(content_gch_HCT116)<-paste(content_gch_HCT116[,1],content_gch_HCT116[,2],content_gch_HCT116[,3],content_gch_HCT116[,4],sep=":")
	
	content_gch_DKO1<-read.table(acc_DKO1,sep="\t",header=F)	
	content_gch_DKO1<-content_gch_DKO1[!duplicated(content_gch_DKO1),]
	rownames(content_gch_DKO1)<-paste(content_gch_DKO1[,1],content_gch_DKO1[,2],content_gch_DKO1[,3],content_gch_DKO1[,4],sep=":")
	
	dataSeq<-seq(5+((length(content_hcg_HCT116[1,])-5)/2)-as.integer(scale/bin_size_align), 5+((length(content_hcg_HCT116[1,])-5)/2)+as.integer(scale/bin_size_align), by=as.integer(move_step/bin_size_align))
	value_hcg_HCT116<-NULL
	value_hcg_DKO1<-NULL
	value_gch_HCT116<-NULL
	value_gch_DKO1<-NULL
	
	for(i in dataSeq){
		
		if(floor(i+move_step/(2*bin_size_align)-1)-ceiling(i-move_step/(2*bin_size_align)) <= 0){
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align))]))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
			value_gch_HCT116<-cbind(value_gch_HCT116,as.numeric(content_gch_HCT116[,ceiling(i-move_step/(2*bin_size_align))]))
			value_gch_DKO1<-cbind(value_gch_DKO1,as.numeric(content_gch_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
		}else{
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(rowMeans(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(rowMeans(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))
			value_gch_HCT116<-cbind(value_gch_HCT116,as.numeric(rowMeans(content_gch_HCT116[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
			value_gch_DKO1<-cbind(value_gch_DKO1,as.numeric(rowMeans(content_gch_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
		}
		
	}
	
	rownames(value_hcg_HCT116)<-rownames(content_hcg_HCT116)
	rownames(value_hcg_DKO1)<-rownames(content_hcg_DKO1)
	rownames(value_gch_HCT116)<-rownames(content_gch_HCT116)
	rownames(value_gch_DKO1)<-rownames(content_gch_DKO1)
	
	
	##clustering:
	##1st step: clustering on DNA methylation combination of 2 samples:
	
	valueTmpMethyNoCluster.1<-value_hcg_DKO1[,40:60]
	valueTmpMethyNoCluster.1<-interpNA(valueTmpMethyNoCluster.1[rowSums(!is.na(valueTmpMethyNoCluster.1))>=3,], method = "linear")
	valueTmpMethyNoCluster.1<-removeNA(valueTmpMethyNoCluster.1)
	
	valueTmpMethyNoCluster.2<-value_hcg_HCT116[,40:60]
	valueTmpMethyNoCluster.2<-interpNA(valueTmpMethyNoCluster.2[rowSums(!is.na(valueTmpMethyNoCluster.2))>=3,], method = "linear")
	valueTmpMethyNoCluster.2<-removeNA(valueTmpMethyNoCluster.2)
	common.methy<-intersect(rownames(valueTmpMethyNoCluster.1),rownames(valueTmpMethyNoCluster.2))
	valueTmpMethyNoCluster<-cbind(valueTmpMethyNoCluster.1[common.methy,],valueTmpMethyNoCluster.2[common.methy,])
	#cluster.methy<-hclust(dist(valueTmpMethyNoCluster),method="ward")
	#order.methy<-rownames(valueTmpMethyNoCluster[cluster.methy$order,])
	##k-means in 2
	cluster.methy<-kmeans(valueTmpMethyNoCluster,methySubClusterNum)
	subClusterOrder<-cluster.methy$cluster
	names(subClusterOrder)<-rownames(valueTmpMethyNoCluster)
	
	##2nd step
	combinedRowOrder<-NULL
	for(clusterNum in c(1:methySubClusterNum)){
		cluster_name<-rownames(valueTmpMethyNoCluster[cluster.methy$cluster == clusterNum,])
		valueTmpCluster.1<-value_gch_DKO1[cluster_name,][,regionToCluster]
		valueTmpCluster.1<-interpNA(valueTmpCluster.1[rowSums(!is.na(valueTmpCluster.1))>=3,], method = "linear")
		valueTmpCluster.1<-removeNA(valueTmpCluster.1)
		valueTmpCluster.2<-value_gch_HCT116[cluster_name,][,regionToCluster]
		valueTmpCluster.2<-interpNA(valueTmpCluster.2[rowSums(!is.na(valueTmpCluster.2))>=3,], method = "linear")
		valueTmpCluster.2<-removeNA(valueTmpCluster.2)
		
		
		common.acc<-intersect(rownames(valueTmpCluster.1),rownames(valueTmpCluster.2))
		valueTmpAccNoCluster<-cbind(valueTmpCluster.1[common.acc,],valueTmpCluster.2[common.acc,])
		cluster.acc<-hclust(dist(valueTmpAccNoCluster),method="ward")
		combinedRowOrder<-c(combinedRowOrder,rownames(valueTmpAccNoCluster[cluster.acc$order,]))
	}
	
	
	returnResult<-list(rowNames=combinedRowOrder, subClusterOrder=subClusterOrder, capUplimit=NULL, capDownLimit=NULL)
	return(returnResult)
	
	
}


clusteringNOMeSeqMatrix2step<-function(methy_HCT116,methy_DKO1,acc_DKO1, accSubClusterNum=4, methySubClusterNum=2){
	content_hcg_HCT116<-read.table(methy_HCT116,sep="\t",header=F)	
	content_hcg_HCT116<-content_hcg_HCT116[!duplicated(content_hcg_HCT116),]
	rownames(content_hcg_HCT116)<-paste(content_hcg_HCT116[,1],content_hcg_HCT116[,2],content_hcg_HCT116[,3],content_hcg_HCT116[,4],sep=":")
	
	content_hcg_DKO1<-read.table(methy_DKO1,sep="\t",header=F)	
	content_hcg_DKO1<-content_hcg_DKO1[!duplicated(content_hcg_DKO1),]
	rownames(content_hcg_DKO1)<-paste(content_hcg_DKO1[,1],content_hcg_DKO1[,2],content_hcg_DKO1[,3],content_hcg_DKO1[,4],sep=":")
	
	content_gch_DKO1<-read.table(acc_DKO1,sep="\t",header=F)	
	content_gch_DKO1<-content_gch_DKO1[!duplicated(content_gch_DKO1),]
	rownames(content_gch_DKO1)<-paste(content_gch_DKO1[,1],content_gch_DKO1[,2],content_gch_DKO1[,3],content_gch_DKO1[,4],sep=":")
	
	dataSeq<-seq(5+((length(content_hcg_HCT116[1,])-5)/2)-as.integer(scale/bin_size_align), 5+((length(content_hcg_HCT116[1,])-5)/2)+as.integer(scale/bin_size_align), by=as.integer(move_step/bin_size_align))
	value_hcg_HCT116<-NULL
	value_hcg_DKO1<-NULL
	value_gch_DKO1<-NULL
	
	for(i in dataSeq){
		
		if(floor(i+move_step/(2*bin_size_align)-1)-ceiling(i-move_step/(2*bin_size_align)) <= 0){
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align))]))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
			value_gch_DKO1<-cbind(value_gch_DKO1,as.numeric(content_gch_DKO1[,ceiling(i-move_step/(2*bin_size_align))]))
		}else{
			value_hcg_HCT116<-cbind(value_hcg_HCT116,as.numeric(rowMeans(content_hcg_HCT116[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
			value_hcg_DKO1<-cbind(value_hcg_DKO1,as.numeric(rowMeans(content_hcg_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))
			value_gch_DKO1<-cbind(value_gch_DKO1,as.numeric(rowMeans(content_gch_DKO1[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
		}
		
	}
	
	rownames(value_hcg_HCT116)<-rownames(content_hcg_HCT116)
	rownames(value_hcg_DKO1)<-rownames(content_hcg_DKO1)
	rownames(value_gch_DKO1)<-rownames(content_gch_DKO1)
	
	
	##clustering:
	##1st step: clustering on DNA methylation combination of 2 samples:
	
	valueTmpMethyNoCluster.1<-value_hcg_DKO1[,40:60]
	valueTmpMethyNoCluster.1<-interpNA(valueTmpMethyNoCluster.1[rowSums(!is.na(valueTmpMethyNoCluster.1))>=3,], method = "linear")
	valueTmpMethyNoCluster.1<-removeNA(valueTmpMethyNoCluster.1)
	
	valueTmpMethyNoCluster.2<-value_hcg_HCT116[,40:60]
	valueTmpMethyNoCluster.2<-interpNA(valueTmpMethyNoCluster.2[rowSums(!is.na(valueTmpMethyNoCluster.2))>=3,], method = "linear")
	valueTmpMethyNoCluster.2<-removeNA(valueTmpMethyNoCluster.2)
	common.methy<-intersect(rownames(valueTmpMethyNoCluster.1),rownames(valueTmpMethyNoCluster.2))
	valueTmpMethyNoCluster<-cbind(valueTmpMethyNoCluster.1[common.methy,],valueTmpMethyNoCluster.2[common.methy,])
	#cluster.methy<-hclust(dist(valueTmpMethyNoCluster),method="ward")
	#order.methy<-rownames(valueTmpMethyNoCluster[cluster.methy$order,])
	##k-means in 2
	cluster.methy<-kmeans(valueTmpMethyNoCluster,methySubClusterNum)
	subClusterOrder<-cluster.methy$cluster
	names(subClusterOrder)<-rownames(valueTmpMethyNoCluster)
	
	##2nd step
	combinedRowOrder<-NULL
	for(clusterNum in c(1:methySubClusterNum)){
		cluster_name<-rownames(valueTmpMethyNoCluster[cluster.methy$cluster == clusterNum,])
		valueTmpCluster.1<-value_gch_DKO1[cluster_name,][,regionToCluster]
		#valueTmpCluster.1<-interpNA(valueTmpCluster.1[rowSums(!is.na(valueTmpCluster.1))>=3,], method = "linear")
		#valueTmpCluster.1<-removeNA(valueTmpCluster.1)
		#maxNAallow=quantile(rowSums(is.na(valueTmpCluster.1)),probs=seq(0,1,0.05),na.rm =T)[19]
		maxNAallow=as.integer(length(regionToCluster)/2)
		valueTmpCluster.1<-valueTmpCluster.1[rowSums(is.na(valueTmpCluster.1))<=maxNAallow,]
		
		valueTmpAccNoCluster<-valueTmpCluster.1
		cluster.acc<-hclust(dist(valueTmpAccNoCluster),method="ward")
		combinedRowOrder<-c(combinedRowOrder,rownames(valueTmpAccNoCluster[cluster.acc$order,]))
	}
	
	
	returnResult<-list(rowNames=combinedRowOrder, subClusterOrder=subClusterOrder, capUplimit=NULL, capDownLimit=NULL)
	return(returnResult)
	
	
}
